From Mott insulator to band insulator: a DMFT study. 
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The question if a Mott insulator and a band insulator are fundamentally different has been the matter of 
intensive research recently. Here we consider a simple model which allows by tuning one parameter to go 
continously from a Mott insulator to band insulator. The model consists of two Hubbard systems connected 
by single particle hopping. The Hubbard Hamiltonian is solved by the Dynamical Mean-Field theory using 
Quantum Monte Carlo to solve the resulting the quantum impurity problem. The quasiparticle spectral function 
is calculated. Here we focus on the optical conductivity and in particular on the Drude weight which can be 
experimentally measured. From our calculation we conclude that there is a continous crossover from the band 
insulator to the Mott insulator phase at finite temperature. 
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I. INTRODUCTION 

Recently the question whether a Mott insulator and 
a band insulator are fundamentally different has been 
raised'&M:Si5ii^. To study this question, we consider the sim- 
plest model which allows, by tuning one parameter, to ob- 
tain a Mott insulator as well as a band insulator phase. The 
model consists of two Hubbard systems with strong on-site 
Coulomb repulsion which are connected by single particle 
hopping. This model can be viewed as a model for two planes 
of strongly correlated electrons on a square lattice with on-site 
interaction and a hopping connecting corresponding sites of 
the two planes. At half filling with no interaction the metal- 
to-band insulator transition is driven by increasing the hop- 
ping between the two subsystems, i. e., the splitting of the 
bonding and antibonding bands produces a gap. In the case of 
no hopping between the planes, a Mott transition is driven by 
increasing the on-site Coulomb repulsion which localizes the 
electrons by suppressing the hopping between different sites. 
The overall scale of the problem is set by the hopping matrix 
element t within the plane (which we set to unity) and the pa- 
rameter differentiating between the Mott and band insulator is 
the ratio of the hopping matrix element between the planes, 
tj^, and the on-site Coulomb repulsion U . The approximation 
used in this work consists of letting the coordination number 
of the sites in each plane (4) to go to infinity. This model has 
been studied with Dynamical Mean-Field theory (DMFT) ap- 
proximation using iterated perturbation theory (IPT) at zero 
temperature''"^. We study this model at finite temperature 



using QMC"''^ as impurity solver Also, first successful at- 
tempts to apply a more demanding continuous-time QMC al- 
gorithm to a simplified two-impurity problem already exist'^. 
Here, the focus is on the nature of the transition from the 
Mott insulator to the band insulator phase. We calculate the 
optical conductivity for direct comparison with experimental 
data. For other theoretical studies of low-dimensional coupled 
strongly correlated systems see e.g. Essler and Tsvelik ' , Pott- 
hoff and Nolting Biermann et al. Koga et al. 

The rest of this paper is organized as follows: First, we in- 
troduce the Hamiltonian and present the solution method us- 
ing DMFT and the Quantum Monte Carlo algorithm. We use 



these methods to determine the metal-to-insulator transition 
and calculate the phase diagram of the model at finite temper- 
ature. Then we consider and analyze the numerical results, 
in particular the analytic continuation of the imaginary-time 
QMC data. We present detailed results for the single-particle 
spectral function and the optical conductivity close to the tran- 
sition and analyze the behavior of these properties close to the 
transition. Finally, we state our conclusions. 



II. FORMALISM 

A. The model and solution method 

The two-plane Hubbard model with interplane hopping t±^ 
is described by the Hamiltonian 



tl-^ctaaC'tc^^l-c (1) 



\,j)aa 



with Ciaa denoting the annihilation operator for an elec- 
tron/hole with spin component a on site i of the plane a = 
0, 1, and rijcTQ = c^o-Cio-a. This means electrons can move 
inside the planes as well as between corresponding sites on 
the two planes, z denotes the coordination number of the lat- 
tice (z = 4 for two dimensions), ensuring a constant band 
width as the coordination number is taken towards infinity. 

Using Dynamical Mean-Field theor y'^i^^'i^' , the two-plane 
system is reduced to two impurities self-consistently embed- 
ded in a bath: In order to calculate on-site (local) properties 
of the sites (i, a) = {i, 0), («, 1) (site i of each of the planes), 
the self-energy ^aa' i^^, k) is replaced by the local self-energy 
Sqq' {lu, 0), leading to a two-impurity (i — 0, a — 0, 1) prob- 
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lem given by the "effective" action 

= J drdr' J2 4.aWeaa'(r,r')-^co.a'(r') 

aa.a' 

+U / dT^noQt(T)rioax(T). (2a) 

The Weiss field Q describes the dynamics of the site i = 
without the interaction plus the rest of the lattice. is a 
2x2 matrix; since the system is symmetric under exchange 
of the planes, we use Gqq = Gn =: Go, Gai = ^lo Gi', 
the properties of the system can be described by the sym- 
metric/antisymmetric (bonding/antibonding) combinations of 
the two planes. This impurity problem is defined by the self- 
consistency equation 



Gs/A {iiOn) 



(2b) 



+i)(iw„ + ^ =F ^± - Ss/A (it^n)) , 



where D{C,) ~ J de D{e) {( — e)^^, D being the density 
of states (DOS) for a free (U = 0) single-plane system, 
Gs/A = ^0 ± Gi, Ss/A = So ± Si, and E is the self-energy for 
the impurity problem, which can be calculated from the ef- 
fective impurity action via the impurity Green's function (the 
mean-field approximation for the on-site lattice Green's func- 
tion). The self-sonsistency equation can be derived exactly 
following the lines given in the work by Georges et al.— . 

Moreover, D is the only place where the detailed lattice 
structure enters the calculations, so the results are essentially 
independent of those details. 

The DMFT equations are usually solved using an iteration 
algorithm consisting of two parts: By solving an impurity-like 
problem ( l2al . the on-site Green's function is determined, then, 
using the DMFT self-consistency equation ( 12b > . a new impu- 
rity problem is defined. This is repeated until convergence has 
apparently been reached. 

We solve the two impurity problem using the Quantum 
Monte Carlo algorithm developed by Hirsch and Fye'^. In or- 
der to use the Monte Carlo algorithm with the DMFT effective 
action which is non-local with respect to imaginary time r, 
the action S has to be rewritten^' using a lattice Hamiltonian 
consisting of auxiliary "bath" orbitals, replacing the "bath" 
Green's function G- 

For initialization, we use a guess for the Weiss fields, 
Gq^^'^'' (ioJn) (diagonal, i. e., connecting one site to itself) and 
Gf'^^iiujn) (off-diagonal, i. e., connecting one site to the cor- 
responding site on the other plane), determining the Green's 
function of the lattice. Using the QMC algorithm, the lo- 
cal imaginary-time Green's functions Go{t) and G'i(t) are 
calculated. Go being the on-site Green's function, whereas 
Gi is again connecting two corresponding sites on differ- 
ent planes. Use of the Dyson equation then yields the self- 
energies Y,o{ii^n) and Ei(ia;„). 

Now, in order to use the self-consistency equation, we 
switch to the symmetric/antisymmetric combinations of the 
two planes, so the self-energy, the Green's function, and the 



Weiss field become diagonal 2x2 matrices. Since the kinetic 
energy is then diagonal as well, the free Green's functions are 
the Hilbert transforms of the density of states for the symmet- 
ric/antisymmetric combinations of the real-space planes with- 
out interaction: 

Glj,{iuJn) = D{tuj„ + ^lTtl_) . (3) 

Therefore the self-energy can be easily calculated using 

Ss/A («W„) = D {iuJri + A* T t±y^ - Gs/A (it^n)"^ , (4) 

where Gs/a = Go ± Gi. Now we can calculate the new Weiss 
fields for the next iteration using the self-consistency equation 

B. Optical conductivity 

Using the electron spectral densities As/\{uj), we calcu- 
late the electron self-energy at real frequencies, J^s/Ai'^)- 
The spectral function for a non-vanishing momentum is 
then given to be A^'^(w) = -lmGs/A{iLJ + i0,e)/7r = 
— Im (1/ + lO =F ij_ — e — Es/Al^i^))) /tt, where e is the 
free-particle kinetic energy. 

The optical conductivity is, up to a constant, defined by 



(5) 



where Gjj denotes the current-current correlation function. 
As a function of the bosonic Matsubara frequencies i/,„, in 
DMFT for a hypercubic lattice, it is— 



G.jj{iiyra) 



E 

Q=S.A 



E 

n— — oo 



deD{e)x (6) 

Ga{iuJn,£)Ga{iuJn + W-m, e)- 



After continuation to real frequencies uj, the real (non- 
dissipative) part of the optical conductivity is, up to the 
constant^^ CTo, 



Re(T(cj) = (To / deD{e) ^ / dcj' 



(7) 



,f{u')-f{u + uj') 



f denoting the Fermi function, f{uj) — 1/ (exp {/Suj) + 1), 
where (3 denotes the inverse temperature. Finally, the weight 
of the Drude peak was determined by fitting a Lorentz curve 
to the central peak (the first five data points, corresponding to 

UJ < 2/(3). 

The reason for the summation over symmetric and antisym- 
metric planes in eqns. (|6j and is the following: since the 
optical conductivity is defined as a long-wavelength limit, the 
momentum transferred by the optical conductivity has to van- 
ish, viz. the in-plane component as well as the component 
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perpendicular to the planes. The perpendicular component 
can assume just two values, and 7r/(plane distance), cor- 
responding to symmetric and antisymmetric orbitals, respec- 
tively. Therefore, the optical conductivity at vanishing (also 
perpendicular) momentum is given by the product of Green's 
functions both symmetric or both antisymmetric. — In the 
limit of high dimension, inter-band transitions do not con- 
tribute to the optical conductivity, as they may only arise from 
interaction vertices. However, in that limit, the interaction 
only contributes to the optical conductivity via self-energy 
insertions in the single-particle Green's function^', quite re- 
gardless of the detailed band-index structure of the interaction 
vertex. 

The replacement of the Gaussian free DOS of the hypercu- 
bic lattice by a semicircular one (which is the exact DOS for 
an infinite-coordination Bethe lattice) is an ad-hoc approxima- 
tion, which may be justified by the low weight of the Gaussian 
tails and their unphysicality. However, eq. (|6} was derived 
for a hypercubic lattice with Gaussian DOS, so, although we 
assume our results to be realistic, they still are based on a sub- 
stantial approximation. - For a detailed discussion of possible 
transport properties on a Bethe lattice, please refer to Bliimer 
and van Donger^^i— . 



in. NUMERICAL RESULTS 
A. Single particle density of states 

We consider the half-filled model (n — 1) at a temperature 
T = 0.025 = 1//3, using L = 100 time slices of At = 
P/L = 0.4. As the density of states of the free {U = 0) 
uncoupled {t± = 0) lattice, we use a semicircular D{£) = 
•\/4 — / (27r), which becomes exact for electrons on a Bethe 
lattice^'. This is more convenient than a Gaussian DOS for 
a hypercubic lattice, because the extended unphysical tails of 
the Gauss distribution render it impossible to clearly define 
the metal-to-band insulator transition. 

From the imaginary-time Green's functions produced by 
the QMC algorithm, the corresponding spectral densities for 
the symmetric/antisymmetric planes are extracted using the 
Maximum-Entropy method^^-^^. We use a default model con- 
sisting of normalized semi-ellipses of half-width C//2 + 2 
centered at =Ft^ for the symmetric and antisymmetric plane, 
resp., plus a small flat "background" in order to keep the pos- 
sibility to extract features outside this area. - Alternatively, a 
flat and a Gaussian default model were used; however, those 
produce unphysically large high-frequency tails in the spec- 
tral density and, as well, artificial humps at = even in 
the non-interacting case. Figure[2shows the density of states 
of the uncoupled system {t^ = 0) at [/ = 1, 2, 3, 4.5, 4.75, 
and 6. The DOS siU — 4.5 has a three-peak shape charac- 
teristic for the metallic state close to the Mott transition. The 
optical conductivity and the Drude weight yield a transition 
value U w 4.7 (Fig. 1^. As the iteration was initialized us- 
ing an "insulating" Green's function, the transition marks the 
lower-U end of the coexistence region. The spectral density at 
f/ = 4.75 represents the insulating state just after the vanish- 




FIG. 1: (Color online) Mott transition. Spectral densities of the un- 
coupled (fx = 0) two-plane Hubbard model. The DOS correspond 
to the metallic state ([/ = 1, 2, 3, 4.5), and to the insulating state 
(U = 4.75, 6). The DOS at [/ = 4.75 corresponds to the insulating 
state slightly above the Mott transition. The iteration was initialized 
using an insulating Green's function. 




FIG. 2: (Color online) Band transition. Reconstructed spectral den- 
sities of the symmetric plane of the free two-plane Hubbard model 
at (7 = 0, = 0, 1,2,3. The antisymmetric DOS is Aa(io) = 
As{—i-o). As and Aa do not overlap in the band insulating state, the 
band transition occurs at t±_ = 2. 



ing of the quasiparticle peak. The DOS at [/ = 6 displays the 
lower and upper Hubbard bands at — J7/2 and +J7/2, resp. 

At [7 = 0, the metal-to-band insulator transition was found 
at tj^ = 2.0 (see also the phase diagram. Fig. 0. This is the 
point where the overlap of the spectral densities for the sym- 
metric and the antisymmetric planes vanishes. The spectral 
densities for the metallic and the band insulating phases are 
given in Fig. |2 As can be seen, the spectral densities for the 
symmetric/antisymmetric plane are shifted from their t = 
position by exactly =Fi^. The symmetric DOS at = 2 cor- 
responds to the state right at the band transition. The error 
bars are of the order of magnitude of the line width. 

For finite U, on increasing t^, for the symmetric plane, the 



4 



' 1 ' 1 ' 1 ' 1 


U = 4.0,t^= 1.4 

\ 

\ 

N 

N. 


— M h-i 

/' P 
/' \ 


U = 4.0, tj^ = 0.4 

\ \ 
\ \ 

, \ 1 ^ , 



0.08 



0.06 




CO 




0.04 



0.02 



FIG. 3: (Color online) Reconstructed symmetric (solid line) and anti- 
symmetric (broken line) spectral densities of the two-plane Hubbard 
model at (7 = 4 and tx_ = 0.4, 1.4 using L = 100 time slices. The 
changes of the Hubbard bands due to tx can be clearly seen. — Due 
to particle-hole symmetry of the two-plane system at half filling, the 
overall spectral density is symmetric. 



weight of the upper Hubbard band is reduced, whereas the 
lower one increases, until, at > 2, the upper Hubbard band 
has completely vanished. For the antisymmetric plane, the 
upper band is increased at the expense of the lower one. For 
intermediate values, this effect can be clearly seen from Fig.|3] 
Our results are compatible with earlier results for the single- 
plane model found by different methods like QMC or IPT^' 
or NRG^^ (the upper boundary of the coexistence region is 
slightly higher in our case, due to very large time slices). As 
discussed below, our results are also compatible to the quan- 
tities calculated for a two-plane model by Moeller et al.i£, 
although those are zero-temperature data. 



B. Optical conductivity 

At first, we consider the optical conductivity of the uncou- 
pled system (<j^ = 0). Using the spectral densities obtained 
by Maximum-Entropy, we found the optical conductivity for 
different values of U (Fig. |4} using eq. 0. Our results are 
compatible with the single -plane data in Pruschke et al. 

The quasiparticle contribution to conduction is given by the 
weight of the Drude peak located at = 0, thus, an insulating 
system has vanishing Drude weight. With increasing interac- 
tion parameter [/, the Drude peak decreases for all values of 
t^. As well, the growth of the incoherent peak at w w [/ is 
clearly visible. In the inset, the Drude weight is shown as a 
function of IJ . Clearly, the system becomes a Mott insulator 
at [/ « 4.7, if the iteration is initialized with an "insulating" 
Green's function. Figure|5]depicts the optical conductivity for 
different at J7 = 2, again consisting of the Drude peak 
of different weights and an "incoherent" part which consists 
of two peaks, one of them located at w [/, the other one, 
present only in the metallic phase, located at w w U jl. How- 
ever, the latter one is usually smeared too strongly to be seen 



FIG. 4: (Color online) The evolution of the optical conductivity as 
a function lo at different U , t± = 0. Inset: formation of the Drude 
weight with increasing U. As can be seen, the Drude weight vanishes 
at (7 ~ 4.7. Because the iteration was initialized using an "insulat- 
ing" Green's function, this is the lower end of the coexistence region. 




FIG. 5: (Color online) The evolution of the optical conductivity as 
a function ui at different t±, (7 = 2. Inset: formation of the Drude 
weight with increasing t±. As can be seen, the Drude weight van- 
ishes at t± ~ 1.8, where a metal-to-insulator transition takes place. 



clearl}i2a, only for low values of tj_, some traces of this peak 
might be recognized. With increasing t±, the Drude peak van- 
ishes at tj^ w 1.8 for U — 2, indicating the transition to a 
predominantly band insulating state. — The transition value 
was found by linear extrapolation of the squared Drude weight 
(see the inset of Fig. |5}. 



C. Phase diagram 

In Fig. |6j the Drude weights for the different values of 
{t± , U) are shown. In order to find the metal-to-insulator tran- 
sitions we used a linear interpolation of the quadratic value of 
the Drude weight, obtaining the phase diagram given in Fig. 
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FIG. 6: Drude weight D at temperature T — 0.025, the iteration was 
initialized with an "insulating" Green's function. The region D 
can thus be identified as the low-U, low-t^ region, as also depicted 
in Fig.0 




FIG. 7: Finite-temperature (T = 0.025) phase diagram. The lines 
are just a guide to the eye. Inset: the evolution of the weight of the 
lower Hubbard band of the symmetric plane on the dotted line. — By 
comparing to Fig. |6| the metallic region is recognized as the region 
with nonvanishing Drude weight. 



The different regions of the phase diagram could be clearly 
located: as expected, there is a metallic state for low U and 
low t±, which is bounded by a metal-to-band insulator tran- 
sition at t± = 2. For high U, the system is in a Mott insu- 
lating state; the metallic and insulating solutions are both lo- 
cally stable within a coexistence region. As discussed below, 
no clear separation between the Mott and the band insulating 
states were found. 

In order to get some impression of the transition between 
the band insulating and the Mott insulating phase, we calcu- 
lated the weight of the lower Hubbard band (LHB) of the sym- 
metric plane which is defined as do; As{iL;), at the points 
given by the dotted line in Fig. Some of the spectral den- 
sity functions can be seen in Fig. |8] The evolution of the LHB 
weight along the dotted line is shown in the inset in Fig. For 
half filling, the weight of the LHB for a purely band insulating 
phase is unity, for a purely Mott insulating phase, it is close to 
0.5. As can be seen the weight of the LHB does neither show 




FIG. 8: (Color online) Selection of reconstructed symmetric spectral 
densities on the dotted line in Fig. at temperature T — 0.025. — 
A purely Mott insulating state is characterized by a spectral density 
of the symmetric or antisymmetric plane, resp., divided by half into a 
lower and an upper Hubbard band, whereas a purely band insulating 
state means the symmetric band is entirely located below tj = 0, and 
the antisymmetric band entirely above. 



some distinct kink nor vanish from a well-defined point. Thus, 
this quantity does not yield any evidence for a phase transition 
between the Mott-Hubbard and the band insulating phase. 

We find a phase diagram which is clearly compatible to the 
zero-temperature phase diagram in Moeller et al.— , keeping 
in mind that J/Moeiier = U/2 and tabMoeiier = t±/2. However, 
some differences ought to be noticed: the coexistence region 
is found at a lower U value due to the finite temperature (see, 
for comparison, the phase diagram in Georges et al.^', where 
the same scale of U is used as by Moeller et alJS); the co- 
existence region has become smaller as well. A coexistence 
region thus clearly exists at a temperature of T = 0.025; in 
contrast, at T = 0.05, no coexistence region was found any 
more. — This behavior suggests that the critical temperature 
of the Mott transition decreases as t± is increasing, see the 
sketch in Fig.|9] 

The other clear difference is the slope of the transition line 
at low interaction U, close to the band transition. We find, for 
low U, the transition line to be at almost constant t±, whereas 
Moeller et al. find a clear dependence on t±. We assume 
this is due to the temperature: as a finite temperature always 
smoothens a metal-to-band insulator transition, a small inter- 
action driving the system to an insulating state can be com- 
pensated by thermal fluctuations. 

The re-entrance behavior seen in the IPT'" cannot be re- 
solved accurately in our calculation. The general shape of the 
phase boundary however suggests that a re-entrance behav- 
ior does not exist at the temperature considered. To resolve 
this issue definitely, lower temperatures have to be consid- 
ered which are inaccessible to the Hirsch-Fye algorithm for 
the large-C/ case. 

The comparative!)*^ high upper bound of the coexistence 
region may be due to the non-negligible Trotter error in this 
region. The insulating solution will be much less sensitive 
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FIG. 9: Three-dimensional piiase diagram of the two-plane Hubbard 
model, composed of the data in Moeller et al.— (T = 0), Georges 
et al.^' (t± = 0), and this work. The dashed line indicates the shape 
of the coexistence region suggested by the data, which are indicated 
by the full lines; thin lines indicate the metal-to-insulator transitions. 

to an increase in the (t/Ar)^ term neglected by the Trotter 
decomposition of the path integral than the metallic solution, 
which is why the upper bound shifts due to the truncation er- 
ror, but not so much the lower bound, which is also at a lower 
U value, so the second-order U term neglected by the time 
discretization is smaller. 



IV. SUMMARY 

In summary, we have calculated the spectral densities, the 
optical conductivities, and the Drude weights of a two-plane 



Hubbard model at low temperature for different values of the 
inter-plane coupling. We have located the different metal-to- 
insulator transitions; however, no clear transition between the 
Mott insulating phase and the band insulating phase could 
be found; as well, the corresponding spectral weights show 
a continuous behavior. This observation is consistent with the 
assumption that there exists only a crossover between those 
two insulating phases, but no clear phase transition. The phase 
diagram is slightly different from the one found in Moeller 
et al. which comes as no surprise as we are considering a 
finite temperature shifting transition values and decreasing the 
coexistence region. 

We have discussed in detail spectral properties like, e. g., 
the optical conductivity, which is of some use to experimental- 
ists. Also, the use of a Quantum Monte Carlo method means a 
serious technical improvement with respect to earlier studies, 
as it is a numerically exact method without the use of uncon- 
trolled approximations. 

Even though the use of a Quantum MC algorithmJ^ means 
a technical advantage, the behavior of the system for very low 
temperatures could not be considered in this work. We plan 
to investigate lower temperatures using a very recently devel- 
oped continuous-time Quantum Monte Carlo algorithmSSiS, 
yielding the phase diagram at much lower temperatures and 
clarifying the evolution towards zero temperature. 
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